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Abstract 

We seek to approximate a composite function /i(x) = g{f{-x)) with a global 
polynomial. The standard approach chooses points x in the domain of / 
and computes /i(x) at each point, which requires an evaluation of / and 
an evaluation of g. We present a Lanczos-based procedure that implicitly 
approximates g with a polynomial of /. By constructing a quadrature rule 
for the density function of /, we can approximate /i(x) using many fewer 
evaluations of g. The savings is particularly dramatic when g is much more 
expensive than / or the dimension of x is large. We demonstrate this pro- 
cedure with two numerical examples: (i) an exponential function composed 
with a rational function and (ii) a Navier-Stokes model of fluid flow with a 
scalar input parameter that depends on multiple physical quantities. 

Keywords: dimension reduction, Lanczos' method, orthogonal 
polynomials, Gaussian quadrature 



1. Introduction &; Motivation 

Many complex multiphysics models employ composite functions, where 
each member function represents a different physical model. For example. 
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consider a simple chemical reaction model; the decay of the concentration 
depends on the decay rate parameter, but the model for the decay rate (i.e., 
the Arrhenius model) depends on the temperature, the gas constant, the 
activation energy, and the prefactor. The concentration is then a function 
of the inputs to the rate model. Studying the effects of these inputs on the 
concentration could be challenging due to the large number of combinations 
of the four rate model inputs. However, one may make use of the composite 
structure - namely, that the concentration depends on these four inputs only 
through the rate model - to find a set of values for the rate parameter with 
which to study the effects on the concentration. By taking advantage of this 
structure, one can hope to reduce the cost of such a study by using fewer 
evaluations of the models. 

We consider the general setting of a composite function 

/i(x) = (7(/(x)), X = (xi, . . . , Xrf) G A- C M'^ (1) 

where 

f -.X — > T C R (2) 

g:J^ — > g C R. (3) 

One may be interested in understanding how h behaves as x changes. How- 
ever, if evaluating h is computationally expensive, then studies that require 
many evaluations may be infeasible. A common approach to this situation is 
to construct a surrogate approximation of h whose evaluations are cheaper 
given inputs x. With a fixed number of h evaluations, one determines the 
parameters (e.g., coefficients) of the surrogate model, and then the surrogate 
is used to study the behavior of h with respect to changes in x. 

Global polynomials [H [2] in x are popular models for surrogate functions 
due to their rapid convergence rate as more terms are added to the approxi- 
mation. For smooth functions, this often translates to relatively few evalua- 
tions to construct an accurate approximation. The approximation properties 
of univariate polynomials {d = 1) are well-studied. Multivariate {d > 1) poly- 
nomial approximations are less well understood, but a simple tensor product 
construction extends the univariate approximation characteristics to the mul- 
tivariate setting. Unfortunately, such tensor product approximations suffer 
from the so-called curse of dimensionality. Loosely speaking, the number 
of function evaluations required to construct an accurate approximation in- 
creases exponentially as the dimension of x increases. When evaluations of 
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h are expensive, this often precludes the use of a multivariate polynomial 
surrogate. 

In this work, we propose a strategy that takes advantage of the composite 
structure of h to reduce the overall cost of building a multivariate polynomial 
surrogate. Suppose one needs m points in X to construct a multivariate 
polynomial surrogate of /i(x), where m may be an exponential function of d. 
Then each evaluation of /i(x) requires first an evaluation of /(x) followed by 
an evaluation of g{f), i.e., 2m total function evaluations. 

Our proposed strategy computes the same m evaluations of /, but uses 
them to implicitly approximate a univariate density function on the range 
space J-'. We then find k m points in J-" on which to evaluate g. The 
process of finding the k points in J-" yields a transformation from the k eval- 
uations of g to m approximations of h; these m approximations of h are then 
used to compute the coefficients of a multivariate polynomial surrogate of 
h. We therefore reduce the cost of the construction from m evaluations of / 
and m evaluations of g to m evaluations of / and k <^ m evaluations of g, 
plus the cost of finding the k points and computing/ applying the transforma- 
tion. This is particularly advantageous when g{f) is much more expensive 
to evaluate than /(x). 

The k points that we find in J-" are the Gaussian quadrature points from 
the implicitly approximated density function of /. The process constructs 
a set of polynomials in / that are orthonormal with respect to the density 
function of /. The function g is then approximated as a truncated series in 
these basis polynomials of /. 

We use a discrete Stieltjes procedure |3] to compute the recurrence coeffi- 
cients of the orthogonal polynomials in /, and we show how this is equivalent 
to a Lanczos' method |3l E] on a diagonal matrix of the evaluations of / with 
a properly chosen starting vector. The basis vectors from the Lanczos itera- 
tion are used to linearly map the k evaluations of g{f) to m approximations 
of g{f{x.)), which are then used to approximate the coefficients of the multi- 
variate polynomial surrogate of h. 

A similar measure transformation approach was used in [SI E] to approx- 
imate functions with sharp gradients. However, the theoretical and compu- 
tational advantages gained by exploiting the connection to Lanczos' method 
were not explored. 

In what follows, we review the polynomial approximation, Gaussian quadra- 
ture, and Lanczos' method (Section [2]); define the problem and derive the 
approximation method (Section^; and demonstrate its applicability on two 
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numerical examples: (i) an exponential function composed with a rational 
function and (ii) a Navier-Stokes model of fluid flow with a scalar input 
parameter that depends on multiple physical quantities (Section |4]). 

2. Preliminaries 

In this section, we briefly review Gaussian quadrature and polynomial 
approximation, as well as Stieltjes' and Lanczos' methods. This will also 
serve to set up notation; we will use notation very similar to Gautschi [3]. 

2.1. Orthogonal polynomials, Gaussian quadrature, and pseudo spectral ap- 
proximation 

Let T C M be equipped with a positive, normalized weight function 
a; : T — 7- IR+ with finite moments; denote an element of T by t. For functions 
u{t) and f (t), define the inner product 

u{t) v{t) uj{t) dt (4) 

r 

with associated norm = (m,'u)-'^/^. Let {7fi(t)} be the set of monic poly- 
nomials (i.e., leading coefficient is 1) that are orthogonal with respect to 
uit), 

^''^'''^■^ = { """o" ' Otherwise. 
The monic orthogonal polynomials satisfy the recurrence relationship 

7r,+i(t) = {t- ai)7Ti{t) - (5iTii-i{t), z = 0, 1, 2, . . . , (6) 

with vr_i(t) = and 7fo(t) = 1. The and /3j are given by 

= —— — , z = 0,1,2,..., (7) 

^ = J^^^^y Z = l,2,.... (8) 

It is often more convenient to work with orthonormal instead of monic or- 
thogonal polynomials, which we write as 

= fl- (9) 
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0,1,2,..., (10) 



The recurrence relationship for the orthonormal polynomials becomes 

with vr_i(t) = and 7ro(t) = 1. If we consider only the first n equations, then 
tTTiit) = \/%ni_i{t) + aiTTi^t) + V'/3i+i7rj+i(t), z = 0, 1, . . . , - 1. (11) 



Setting 7r(t) = [7ro(t), 7ri(t), . . . ,7r„_i(t)]^, we can write this conveniently in 
matrix form as 

t7r(t) = J7r(t) + v^7r„(t)e„, (12) 

where e„ is a vector of zeros with a one in the last entry, and J = (known 
as the Jacobi matrix) is an n x symmetric, tridiagonal matrix 



(13) 



The zeros of iTnit) are the eigenvalues of J and 7r(Aj) are the correspond- 
ing eigenvectors; this follows directly from (12). Let Q be the orthogonal 
matrix of eigenvectors of J; the elements of Q are given by 



k(A,)ll2 



z,j = 0, . . . - 1, 



(14) 



where || ■ II2 is the standard 2-norm on M". We write the eigenvalue decom- 
position of J as 

J = QAQ^. (15) 

It is known that the eigenvalues {Xj} are the nodes of the n-point Gaussian 
quadrature rule associated with the weight function w(t). The quadrature 
weight Uj corresponding to is equal to the square of the first component 
of the eigenvector associated with A-,-, 
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lk(A,)lli' 



(16) 



The weights {uj} are known to be strictly positive. It will be notationally 
convenient to define the matrix W = diag([y^, . . . , ^z/.„_i]). 
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For an integrable scalar function f{t), we can approximate its integral 
by an ra-point Gaussian quadrature rule. Let fj = f{Xj). The quadrature 
approximation is a weighted sum of function evaluations, 

« n—l 

/ f{t)oo{t)dt = J2f,'^j + Rif)- (17) 

■^'^ 3=0 

If f{t) is a polynomial of degree less than or equal to 2n — 1, then R{f) = 0; 
that is to say the degree of exactness of the Gaussian quadrature rule is 
2n-l. 

A square integrable function f{t) admits a mean-squared convergent Fourier 
series in the orthonormal polynomials. A pseudospectral approximation of 
f(t) is constructed by first truncating its Fourier series at n terms and ap- 
proximating each Fourier coefficient with a quadrature rule. If we use the 
n-point Gaussian quadrature, then we can write 

n-1 



fit) ^ 5^/.7r,(t), (18) 

where 



i=0 



n-1 



/, = 5^/,.7r,(A,)z.,. (19) 

j=0 



Let f = [/o, . . . , fn-iV and f = [/o, . . . , fn-iV- Then using Q from (fTsl) 



(18), and (19), we can write 

f = QWf fit) ^ F7r(t) = f'^W^Q'^7r(t). (20) 

The expression f = QWf is the discrete Fourier transform for the orthogo- 
nal polynomial basis. Note that it is easy to show that the pseudospectral 
approximation interpolates fit) at the Gaussian quadrature points. 

2.1.1. Tensor product extensions 

The above concepts extend to multivariate functions via a tensor product 
construction. Let T = 7i ® • ■ • ® 7^ be the domain with elements t = 
(ti, . . . , trf). We assume that the weight function w : T — t- IR+ is separable, i.e. 
a;(t) = oJiitt) ■ ■ -ooditd), where each univariate weight function is normalized 
and positive with finite moments. 
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Tensor product Gaussian quadrature rules are constructed by taking cross 
products of univariate Gaussian quadrature rules: 

\l,...,id = (-^ii' • • • ' -^id)' (21) 

where the points Aj^ with ir = 0, . . . , — 1 are the univariate quadrature 
points for Uritr), with r = 1, . . . ,d. The associated quadrature weights are 
given by the products 

= i^h ■ ■ ■ i^i,- (22) 
To approximate the integral of fit), compute 

/ni-l "d-l 
Jit) uit) dt ^ /(^n,.... J ^n,.....- (23) 

n=0 id=0 

The tensor product pseudospectral approximation is given by 

"1-1 "d-i 

/(t) ^ E---E/^-.^<^^'^(^i)---^^<*(^'^)' (24) 

il=0 id=0 

where 

ni-l rid-l 

fh,...,ia = XI ■ ■ ■ XI /('^n,-..,»d) T^iii\i) ■ ■ ■ T^iS^d) (25) 

ii=o jd=o 

The matrix notation extends via Kronecker products. Let TZn^itr) be the 
vector of univariate polynomials that are orthonormal with respect to Uritr) 
for r = 1, . . . ,d. Then the vector 

7r(t) = 7r„,(ti)®---®7r„^(td) (26) 

contains multivariate polynomials that are orthonormal with respect to a;(t). 
Similarly, define the matrices 

Q = Qni®---®Qn„ W = ® ■■■®W„,. (27) 

Then 

f = QWf fit) ^ f^7r(t) = f^W^Q^7r(t), (28) 

where f is an (rii ■ ■ ■ naj-vectoi of the tensor product pseudospectral coeffi- 
cients, f is an (ni ■ ■ ■ nrf)-vector containing the evaluations of / at the points 
given by (21), ordered appropriately. Again, the expression f = QWf is the 
(i-dimensional discrete Fourier transform for the polynomial basis. 



7 



2. 2. Stieltjes ' procedure 

Stieltjes proposed a procedure for iteratively constructing a sequence of 
univariate polynomials that are orthogonal with respect to a given measure; 
see p]. His method exploits the recurrence relationship for the orthogonal 
polynomials. He observed that with the weighted inner product Q, one may 
begin with 7r_i and ttq, compute ao and /3i from ^ and ([s]), construct vri 
from ([6]), compute ai and (32, construct tc2, and so on. 

A normalized version of Stieltjes' method for computing the orthonormal 



The 



polynomials and their recurrence coefficients is given in Algorithm 
computed ai from Algorithm [l] are equivalent to the expression in ( 7 ) , and 
the computed rji are equal to in 



Algorithm 1 A Stieltjes procedure for computing the first n orthonormal 
polynomials given a normalized weight function uj{t). Let 7r_i = and 

ttq = 1. 

for i = to n — 1 do 

Vi = llTTill 
TTj = TTi/r]i 

ai = (tni, TTj) 

TTj+l = (t- ai)'Ki - r/jTTi.i 

end for 



Gautschi [3] proposed to use a discrete inner product - e.g., based on a 
Gaussian quadrature rule. In the univariate case [d = 1), this becomes 

m— 1 

{u,v) ^ ^M(A,)t;(A,)z/„ (29) 
i=o 

where A-,- and Vj are the points and weights of the discrete inner product. He 
reasoned that if the discrete inner product converges to the continuous, then 
the recurrence coefficients computed with the discrete inner product will also 
converge. Similarly, one may think of the recurrence coefficients computed 
with the discrete inner product as approximations of those computed with 
the continuous inner product. 

In Section [3| we will use a tensor product quadrature rule to define a 
discrete inner product on the space X. We will use this inner product in a 
Stieltjes procedure to compute the recurrence coefficients of a set of univariate 
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orthogonal polynomials, where orthogonality is with respect to the density 
function of / defined on J-". We have chosen to focus the presentation by 
using the tensor product quadrature rule to define the discrete inner product. 
However, the construction can be easily adjusted to employ other discrete 
inner products. 

2.3. Lanczos method 

Lanczos' method [1] for symmetric matrices is the foundation for itera- 
tive eigensolvers and Krylov subspace methods for solving symmetric linear 
systems. It generates a symmetric, tridiagonal matrix (the Jacobi matrix) 
and a sequence of mutually orthogonal (in exact arithmetic) vectors known 
as the Lanczos vectors. The eigenvalues of the tridiagonal matrix - known 
as the Ritz values - approximate the eigenvalues of the symmetric matrix. 

In fact. Algorithm [l] is exactly a form of Lanczos' methocQ if we replace 
(i) the variable t by a symmetric matrix A of size mxm, (ii) the polynomials 
7f i{t) by the Lanczos vectors Vj, (iii) the starting polynomial ttq by a starting 
vector vo, and (iv) the inner product by a discrete, weighted inner product. 

Suppose that k iterations of the method have been executed. We can 
write the recurrence relationship for the Lanczos vectors in matrix notation 
as 

AV = VT + %Vfce^, (30) 

where V = [vq, . . . , Vfc_i] is an m x matrix of Lanczos vectors, T is the 
k X k symmetric, tridiagonal Jacobi matrix of recurrence coefficients, and 
is a last column of the k x k identity matrix. 

3. The approximation method 

In this section, we take advantage of the relationship between an approx- 
imate Stietjes' procedure with a discrete inner product and Lanczos' method 
to devise a computational method for approximating composite functions. 
Recall the problem setup from the introduction: 

h{x.) = g{f{x.)), x = {xi,...,Xd) ex CR" (31) 



^However, Algorithm 1 has undesirable numerical properties as an implementation. 
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where 



f -.X — > C R (32) 
g-.T — > g C R. (33) 

We restrict our attention to input spaces that are d-dimensional hypercubes, 
X = [-1, 1]"', although this can be relaxed to more general hyperrectangles. 
We assume that X is equipped with a positive, separable weight function 
i^x = (^(x) = coi{xi) ■ ■ -codixd) with finite moments. We also assume that / 
is bounded, so that J-" is a closed interval in M. 

A tensor product pseudospectral approximation of h on the space X fol- 
lows the construction in Section I2.1.1I 

ni— 1 na—l 

^(x) ~ ^■■■'^hiu...M^ii{xi)---T^iA^d), (34) 

h=0 id=0 

where the 7rj^(xr) are the univariate polynomials that are orthogonal with 
respect to LOr{xr)- The coefficients hi^^,,,^i^ are 

ni— 1 rid— I 

K,...,id = J2'"Y1 ^^K,-,u) ■ ■ ■ ^iA^a) ^h,-,u^ (35) 

ii=o id=o 

where \ji,...,jd = • • • 5 -^jd) Vj^,--- ,jd — ^h'''^3d nodes and 

weights, respectively, of the tensor product Gaussian quadrature for the 
weight function u^. This can be written conveniently in matrix form as 



in (|28]) as 

h = QxW.h, /i(x) ^ h^7r(x) = h^W^Q^7r(x), (36) 

where h is a vector of the tensor product pseudospectral coefficients; h is 
a vector of evaluations of h at the quadrature points; 7r(x) is a vector of 
the product type multivariate orthonormal polynomials; and Wx and Qx 
represent the c?-dimensional discrete Fourier transform for functions defined 
on X . 

We assume that the primary expense of this computation comes from the 
evaluations of h at the quadrature points. There are m = ni ■ ■ ■ points 
in the tensor product quadrature rule; if rii = ■■■ = Ud = n, then m = n'^. 
Each evaluation of h requires an evaluation of / followed by an evaluation 
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of g. Our goal is to approximate the values of h at the quadrature points 
using m evaluations of /, but only k <^ m evaluations of g\ we accomplish 
this goal by taking advantage of the composite structure of h. 

Our strategy is to apply Stieltjes' procedure to implicitly construct a set 
of polynomials that are orthonormal with respect to the density function of 
/ defined on if. The /c-point Gaussian quadrature rule for the univariate 
density function provides a set of points in J-" on which to evaluate g. We 
then linearly map the k evaluations of g to approximate the m evaluations 
of h. 

Let Uf : T ^ R+ be the normahzed density function of /. We assumed / 
was bounded on X, which implies J-" is a closed interval and LOf is bounded. 
Therefore, ojf has finite moments, and it admits a set of univariate orthonor- 
mal polynomials 0j = with i = 0, 1, 2, . . . . Then we can approximate 
g = g{f ) with a univariate pseudospectral approximation, 

k-1 

9{f) ^ Y.9^M^ (37) 

i=0 

where 

k-1 

i=o 

The 9j G J-" and iij G M+ are the nodes and weights, respectively, of the 
/c-point Gaussian quadrature rule for a;/. In the matrix notation, 

g = Q/W/g, g{f) ^ g^0(/) = g^WjQj0(/), (39) 

where g is a vector of the pseudospectral coefficients; g is a vector of eval- 
uations of g at the quadrature points Oj; <p{f ) is a vector of the univariate 
orthonormal polynomials; and Wj and Q/ represent the discrete Fourier 
transform for functions defined on J-'. 

Notice that if we were given (i) the coefficients gi and (ii) the evaluations 
of the polynomials 0j at the point /(x) with x G A", then 

fc-i 

/i(x) = g{f{K)) ^ J^diMm)- (40) 

i=0 

Unfortunately, the direct evaluation of this expression is problematic, since 
each new point x requires the computation of /(x). We would prefer a 
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surrogate like (34) whose evaluation is independent of both g and /. With 
this in mind, we seek to evaluate (40) at only the quadrature points Aj^,...,i^ 
on X: 



(41) 



1=0 



We denote the output of the pseudospectral expansion (41) by h. We write 
(41) in matrix notation as 



Ug. 



(42) 



The elements of the m x k matrix U are the evaluations of 0j(/(Ajj^...^j^)), 
where each row of U corresponds to a quadrature point Aij,...,j^ and each col- 
umn corresponds to a polynomial 0j(/). This matrix enables a construction 
similar to (34). 

Putting these ideas together, we construct a global polynomial surrogate 
for h{x.) that approximates the tensor product pseudospectral approximation. 
The form of the polynomial surrogate is the same as the true pseudospectral 
approximation - i.e., a linear combination of the multivariate orthogonal 
basis polynomials - but the pseudospectral coefficients are approximated 
using <C m evaluations of g plus the cost of computing the polynomial 
evaluations 4>i{f{\ii,...,iJ)- We will denote the approximate coefficients by 
h^h, which is consistent with the notation h^h. 

More precisely, we construct a surrogate for /i(x) as 

ni-l rtd-l 

^(x) ~ ^■■■^k,,...,i^T^h{^i)---T^i^{xd), (43) 



ii=0 



id=0 



where 



,3d- 



(44) 



ni-l "d-l 

^ ■ ■ ■ ^ H^h, -Jd) T^hi^ji) ■ ■ ■ ^idi^jd) 

31=0 jd=0 

The quantities /i(Ajj,...j-^) are given by (41). In matrix notation, we combine 
(g, (g, and (|42|) to get 

h{x) ^ h^7r(x) 

= h^W^Q^7r(x) 
= g^U^W^Q^7r(x) 



g^WjQ^U^W^Q^7r(x), 



(45) 
(46) 
(47) 
(48) 
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where h is a vector of the coefficients h, , and h is a vector of the eval- 
nations of h. Notice that g contains k evaluations of g. In what follows, 
we will show how to compute the points 6j in ( p8| and the transformations 
U, Q/, and Wj from (48) using m evaluations f{^ji,...,jd)- Thus, we will 
approximate the tensor product pseudospectral approximation (34) using k 
evaluations of g and m evaluations of /. 

3. 1 . Lanczos ' method for approximation 

The strategy is implemented computationally with Lanczos' method ap- 
plied to a diagonal matrix whose nonzero elements are the m evaluations of 
/ at the quadrature points on X; the A;-point Gaussian quadrature rule on 



J-" comes from the computed Jacobi matrix (see (15)), and the map from k 
evaluations of g to m evaluations of h comes from the Lanczos vectors. 

It is notationally convenient in this section to order the d-dimensional 
quadrature points G A" so as to be indexed by the natural numbers. 

For the remainder, we will refer to a node as Aj with corresponding weight z/j 
for i = 0, . . . ,m—l. The specific ordering must be consistent with the tensor 
product structure of (27). We will similarly order the function evaluations 

^ = f{K). 

We first show that Lanczos' method yields the quantities desired from the 
Stieltjes procedure. The fact has been observed elsewhere in the literature [SI 
El [ini [H] ; we state it as a theorem for reference and notation. 



Theorem 1. Let 



fo 



m—l 



(49) 



be the diagonal matrix whose nonzero elements are the evaluations of f at the 
quadrature nodes Xi E X . For an m-vector of ones e and a starting vector 
vq = Wxe, Lanczos method applied to A is equivalent to Stieltjes' procedure 
with a discrete inner product to construct the recurrence coefficients of the 
polynomials (f)i{f) that are orthonormal with respect to an approximation of 
the measure Uf. The discrete inner product is defined by the nodes Aj and 
weights Vi of the tensor product Gaussian quadrature rule for the measure 

Proof. To prove this statement, we simply describe the quantities in Algo- 
rithm hi with the discrete inner product. The starting polynomial 0o = 1 
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corresponds to an m-vector of ones, e. Let v_i be an m-vector of zeros. 
Then the quantities from Algorithm [T] with the discrete inner product be- 
come 



~T~ \ 1/2 




Vi 



T 



ai 



m— 1 

i=o 

<Pi+i{fj) = {fj - ai)(l)^{fj) ~ r]i(f)i.i{fj), 
which can be written 

Vi+i = (A - ail)vj - ?7iVi_i, (50) 
with Vi+i = [0i+i(/o), . . . To recover the polynomials, 

U = W^^V, (51) 
where U(j,i) = 0i(/j). □ 



The matrix U from (51) is the same as the one from (42). Note that 
- as mentioned in Theorem [T] - these quantities do not correspond to the 
exact measure Uf. They instead correspond to an approximation of from 
the evaluations fj = /(Aj). We will not examine the error made in this 
approximation; we assume the evaluations of / at the nodes Xj G X are 
sufficient to resolve the salient features oi Uf. 

Running k steps of the Lanczos process yields the recurrence relationship 
( 30 ) . The elements of the kxk tridiagonal Jacobi matrix T are the recurrence 
coefficients of the polynomials 0i(/) up to order k — 1. The Lanczos vectors 



V contain the polynomial evaluations (j)i{fj) scaled by ^^zTJ as in (51). 
Denote the eigendecomposition of T by 

T = Q^e Qj, e = diag([0o, . . . , Ok-i]). (52) 

The Ritz values (the eigenvalues of T) are the Gaussian quadrature nodes 
for the approximation of on the space J-", and the weights come from the 
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first component of tlie eigenvectors of T as in (16). Precisely speaking, the 
4>i{f) are orthogonal with respect to the discrete measure defined by the 6j 
and fij. 



Therefore, we compute the quantities 9j, fij, Q/, Wj, and U from (48) 
using k steps of Lanczos' method applied to the diagonal matrix A followed 
by the eigendecomposition of the Jacobi matrix T. This is the computational 
cost incurred beyond the m evaluations of / and k evaluations of g. 

3.2. Loss of orthogonality and stopping criteria 

We have stated that we expect k <^ m, ot that the number of points in 
the discrete measure on J-" will be much smaller than the number of points 
in the discrete measure on X. The number k is the number of iterations of 
the Lanczos procedure; how do we know how many iterations to use to get 
an accurate approximation of the measure on J-"? 

It is well known that Lanczos' method in finite precision behaves differ- 
ently than the algorithm in exact arithmetic; a thorough treatment of this 
subject can be found in Meurant's excellent monograph [5], as well as [T2] . 
In particular, the Lanczos vectors will often lose orthogonality after some 
number of iterations. 

Thanks to the work of Paige [13] as described in [H] - as well as [HI [TSl [IS] 

- we know that the loss of orthogonality is closely related to the convergence 
of the Ritz values to the true eigenvalues; loosely speaking, once a Ritz value 
has converged to an eigenvalue, the remaining Lanczos vectors lose orthog- 
onality. It has been observed that in many cases the extremal Ritz values 
converge to the extremal eigenvalues fastest depending on the starting vector. 
From this we can expect that the Lanczos vectors will lose orthogonality once 
the extremal Ritz values are sufficiently close to the extremal eigenvalues. We 
use this expectation to motivate a heuristic for stopping the Lanczos itera- 
tion. Further justification of the following heuristic is the subject of on-going 
work. 

Recall that A is diagonal, so we are not concerned with any particular 
eigenvalue. In fact, we are only concerned with approximating the range of 
the data - which is the range of the function /(x) evaluated at the points Aj 

- and its corresponding measure. Therefore, once the extremal Ritz values 
converge, we are satisfied. Leveraging the work on Lanczos' method in finite 
precision, we can judge when the extremal Ritz values have converged by 
checking orthogonality of the Lanczos vectors. Essentially, we can treat the 



15 



loss of orthogonality in the Lanczos vectors as stopping criteria. We use the 
following measure of loss of orthogonality given a tolerance TOL: 



r = logio(III-V^Vll^) > TOL (53) 

where || ■ ||f is the Frobenius norm. Other efficient measures for determining 
loss of orthogonality are discussed in [17, Chapter 9] as well as [18], [19]. In 
the numerical examples of Section |4| we choose T0L=-14. 

If the iterations continue beyond this point, we find that the points and 
weights of the quadrature rule for the measure on T become less smooth; 
this phenomenon is similar to choosing the wrong bin size for a histogram. 
In some cases, we observe the familiar (to those who have studied Lanczos' 
method) appearance of ghost eigenvalues. If we examine the weights corre- 
sponding to pairs of nearly identical Ritz values, we usually find that one 
of the weights is orders of magnitude smaller than the other. Of course, we 
would prefer to ignore points with very small weights, since this would cor- 
respond to a wasted function evaluation in the quadrature approximations. 
We demonstrate these phenomena on the following numerical examples. 

3.3. An Algorithm 

We close this section with a step-by-step algorithm using the linear alge- 
bra notation to summarize the procedure. Given functions /(x) and (?(/), 
the goal is to approximate the coefficients of a tensor product pseudospectral 
polynomial surrogate for the composite function /i(x) = g{f(x.)). 

1. Obtain the m nodes Aj and weights Ui of the tensor product Gaussian 
quadrature rule for the space X. Also, obtain the m x m matrix Qx 
and the diagonal matrix of the square root of the quadrature weights 
Wx from 



For i = 0, . . . , m — 1, compute = /(Aj), and form the diagonal matrix 
A = diag([/o,...,/m_i]). 

For the starting vector vq = Wxe, run Lanczos' method until r > TOL 



from (53), and let k be the number of iterations. Store the matrix U 
from (51 ). 



■^In a real implementation, Qx and Wx do not need to be formed explicitly. We only 
need the action of Qx on a vector, which can be computed efBciently using methods such 
as [20]. 
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4. Compute the eigenvalue decomposition of the k x k Jacobi matrix T 
from the Lanczos procedure to get the quadrature nodes {9k} and the 



discrete Fourier transform matrices Q/ and W/; see (52). 

5. For i = 0, . . . , k — 1, compute gi = g{6i) and form the vector g = 

6. Compute approximate coefficients for the pseudospectral approxima- 
tion h h as 

h = QxW,UQ/W/g. (54) 

These coefficients define a polynomial approximation of /i(x) with a basis of 
multivariate product-type orthonormal polynomials. 

4. Numerical Examples 

We present two numerical studies demonstrating the qualities of the 
method. The first is an example with functions chosen to stress the method's 
properties. The second applies the method - as a proof of concept - to a 
model from fiuid dynamics with a scalar input parameter that depends on 
multiple physical quantities. 

4-1. Simple functions 

Let x = {xi,X2) G [—1, 1]^ with a uniform measure of 1/4 in [—1, 1]^ and 
zero otherwise. Given parameters 5i> 1 and ^2 > 1, define the function 

/W = 7 A} FT- (55) 

Notice that /(x) > 0, and 8\ and 82 determine how quickly / grows near 
the boundary. The closer b\ and 62 are to 1, the closer the singularity in the 
function gets to the domain, which determines how large / is at the point 
[x\ = 1,0:2 = 1). For the numerical experiments, we choose 6\ = ^2 = 1-3. 
The function / is analytic in x, so we expect polynomial approximations to 
converge exponentially as the degree of approximation increases. 
Next we choose g{j^ = exp(/), so that 



Again, g{f) is analytic in /, so /i(x) is analytic in x, as well. 
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We choose the discrete measure on A' to be a tensor product Gauss- 
Legendre quadrature rule on [—1,1]^ with n points in each variable, which 
results in m = points and weights. The m x m diagonal matrix A has 
diagonal elements equal to / evaluated at the points of the discrete measure. 

To test the approximation properties of the Lanczos-enabled method, it 
is sufficient to examine the error at the tensor product quadrature points; 



see (41). This is essentially the same as computing the difference between 
the true pseudospectral coefficients h and their approximation h. For the 
number m of bivariate nodes of the tensor product quadrature rule and the 
number k of Lanczos iterations, define the error S = S{m, k) as 

m—l k—1 

= Y,{h{\) - Y,9^cl>,{f{X,))f (57) 

= l|h-Ug||2. (58) 

In Figure [l| we plot both log^Q^S) and the measure of orthogonality of the 
Lanczos vectors (see ([53])) as m and k increase. To read these plots, choose m 
from the y-axis, and then follow the plot to the right to increase the Lanczos 
iteration k. It is interesting to note that the error in the approximation 
does not increase after the Lanczos vectors lose orthogonality. The loss of 
orthogonality is useful for a stopping criteria to determine the smallest k 
that produces an accurate approximation. However, taking more than the 
minimum number of Lanczos iterations does no harm for this example. 

In Figure [2} we plot a series of bar graphs of the quadrature weights ni at 
points 9i for the measure on J-" computed with m = 81. While the bar plot 
resembles a histogram, the comparison between a histogram and quadrature 
weights is not precise. Nevertheless, the series of bar plots demonstrates 
the behavior of the weights as the Lanczos iteration index continues beyond 
the point when the Lanczos vectors lose orthogonality; the orthogonality 



measures from (53) are presented in each plot. We observe that the weights 
lose smoothness as the Lanczos vectors lose orthogonality; note the weights 
in the right tail of the plot. 

Fluid flow example 
As an example of applying these techniques to an engineering problem, 
we examine a simple channel flow problem with a scalar input parameter (the 
Reynolds number) that depends on multiple physical quantities (density and 
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Figure 1: Figure la plots the error in approximation of h as measured by 
(57). Figure lb shows the loss of orthogonality in the Lanczos vectors using 
T from (53). 



viscosity). Consider the two-dimensional rectangular domain of length L = 1 
m and width = 0.1 m shown in Figure |3| Water flows into the left side 
of the domain with a horizontal velocity of Uq = 0.01 m/sec, and we are 
interested in computing the velocity of the flow out of the domain on the 
right side. 

At room temperature and standard pressure, the dynamics of the fluid 
within the domain are well-modeled by the incompressible Navier-Stokes 
equations 

p— + p(u • Vu) - /iV'u + VP = 0, (59) 

P(V ■ u) = 0, (60) 

where u is the two-component velocity of the fluid, p is the density, p is the 
viscosity, and P is the pressure. Using the inlet flow velocity and the width 
W of the domain, the equations are non-dimensionalized resulting in 

^ + u ■ Vu - — V'u + VP = 0, (61) 
ot Re 

V ■ u = 0, (62) 

where V = V /W, and 

Re=P^ (63) 
A* 
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B 0,25 

^ 9 



T^--14.51 
T^=-15.21 




Quadrature Points {e.j 



(a) k 




Quadrature Points 

(b) fc = 10 




x^--7.93 
T^=-9.91 




Quadrature Points {9|} 

(c) fc = 15 



Quadrature Points {9|} 

(d) fc = 20 



Figure 2: A series of bar plots showing the weights fii at the points 6i for 
the quadrature rule on J-". The number r in each plot shows the measure of 



orthogonality in the Lanczos vectors; it is defined in (53) 
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u = 



u = (uo,0) 



L W 
u = 



Figure 3: Fluid flow domain 
is the Reynolds number. 



Equations [6IH62] are discretized spatially on a mesh of 500 by 50 quadran- 
gle cells using the finite element method with piecewise bilinear basis func- 
tions for both the velocities and pressures [21J. Given a Reynolds number, 
the resulting nonlinear algebraic equations are solved via Newton's method 
using a GMRES linear solver [22j and incomplete-LU factorization precon- 
ditioner. The resulting flow solution at density p = po = 998.205 km/m^ 
and viscosity p, = Hq = 0.001001 Ns/m^ is shown in Figure [Zj the density 
and viscosity values roughly correspond to water at room temperature and 
standard pressure. The calculations were implemented in the Albany [23] 
simulation package using numerous solver and discretization packages from 
the Trilinos framework [21]. 

We consider a problem where the ambient temperature and pressure are 
uncertain resulting in uncertain density and viscosity. In particular we model 
the density and viscosity as uniformly distributed random variables 

p e [0.99po, l.Olpo] P e [0.9po, 1.1/io]. (64) 

In other words, we assume density varies uniformly by 1% and viscosity varies 
uniformly by 10%. In the notation of Section |3| we have 

X = (p,/i), 

f( ) = — = ^ 

Re puqW 
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Figure 4: Horizontal flow velocity at mean density and viscosity 

The function h{x.) = g{f(x.)) corresponds to the maximum outflow velocity at 
the right side of the domain given fixed values for p and /i. Each evaluation of 



g involves an expensive solution of equations |6Tp2] - compared to computing 
/(x). 

For this experiment, we choose a tensor product Gauss-Legendre quadra- 
ture rule with 11 points in the range of p and 11 points in the range of p for 
a total of m = 121 points. We use the procedure from Section [3] to approxi- 
mate the maximum outflow velocity at all 121 pairs of {p, p) by constructing 
a /c-point Gaussian quadrature rule for 1/Re with k = 13. In other words, 
with only 13 evaluations of g - the expensive flow solver - we can approx- 
imate the output h at 121 points in the parameter space corresponding to 

X. 

To check the error in the approximation, we also compute the maximum 
outflow velocity at all 121 combinations of p and p, which enables the com- 



putation of (57). With 13 steps of the Lanczos procedure, we have a loss of 



orthogonality in the basis vectors of r = —13.14 (see equation (53)). The 
error in approximation (equation (57)) is 1.55 x 10~^. 



5. Conclusion 

We have presented a method for approximating a composite function by 
implicitly approximating the outer function as a polynomial of the output 
of the inner function. This measure transformation is based on Stieltjes' 
method for generating orthogonal polynomials given an inner product, and 
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it is implemented as Lanczos' method on a diagonal matrix of inner function 
evaluations at the points of a discrete measure. We have developed a heuristic 
for when to terminate the Lanczos iteration based on the loss of orthogonality 
in the Lanczos vectors - a common phenomenon for the algorithm in finite 
precision. The resulting method reduces the number of evaluations of the 
outer function, which are only required at the Gaussian quadrature points of 
the transformed measure. The numerical experiments show the behavior of 
the method and the scale of the reduction. 
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